Author

Scott Forrest and Kyana Pike

Published

March 13, 2025

Here we are plotting the results of integrated step selection function (iSSF) models fitted to water buffalo in Arnhem Land. These models contain the following continuous predictors:

Load packages

Code 
library(tidyverse)
library(glmmTMB)
library(terra)
library(RColorBrewer)
library(lemon)

Import data

Code 
ssf_data <- read_csv("pheno_start_end_dat_NBR_lag_NAFImask_ssf10rs_2024-03-21.csv")
Rows: 1062908 Columns: 58
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr   (4): pheno_start, pheno_end, season, day_period
dbl  (51): ...1, id, burst_, x1_, x2_, y1_, y2_, sl_, ta_, dt_, step_id_, sr...
lgl   (1): case_
dttm  (2): t1_, t2_

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Get scaling factors for the NDVI and NDVI^2 covariates

Code 
ssf_data <- ssf_data %>% mutate(srtm_start_scaled = scale(srtm_start),
                            srtm_end_scaled = scale(srtm_end),
                            water_dist_start_scaled = scale(water_dist_start),
                            water_dist_end_scaled = scale(water_dist_end),
                            ndvi_start_scaled = scale(ndvi_start),
                            ndvi_end_scaled = scale(ndvi_end),
                            nbr_start_scaled = scale(nbr_start),
                            nbr_end_scaled = scale(nbr_end),
                            sl_scaled = scale(sl_),
                            log_sl_scaled = scale(log_sl_),
                            cos_ta_scaled = scale(cos_ta_)
)

# get scaling factors
sl_scale_sd <- attr(ssf_data$sl_scaled, "scaled:scale")
log_sl_scale_sd <- attr(ssf_data$log_sl_scaled, "scaled:scale")
cos_ta_scale_sd <- attr(ssf_data$cos_ta_scaled, "scaled:scale")
ndvi_scale_sd <- attr(ssf_data$ndvi_start_scaled, "scaled:scale")
ndvi_scale_sd <- 0.1262109 # from line of code above

Check the NBR layers and NAFI mask

Import layers

Code 
nbr <- rast("nbr_raster_amt_cropped_epsg32753.tif")
terra::time(nbr) <- as.POSIXct(lubridate::ymd("2018-01-01") + months(0:23))

# the NAFI layers - monthly fires indicated by the value
NAFI2018_cropped_resampled <- rast("NAFI2018_cropped_resampled_raster_epsg32753.tif")
NAFI2019_cropped_resampled <- rast("NAFI2019_cropped_resampled_raster_epsg32753.tif")

nbr2018_NAFImasked <- rast("nbr2018_NAFImask_epsg32753.tif")
nbr2019_NAFImasked <- rast("nbr2019_NAFImask_epsg32753.tif")
nbr_NAFImask <- c(nbr2018_NAFImasked, nbr2019_NAFImasked)

Plot the NAFI layers

Code 
plot(NAFI2018_cropped_resampled, col = rev(brewer.pal(10, "RdBu")))

Code 
plot(NAFI2019_cropped_resampled, col = rev(brewer.pal(10, "RdBu")))

Plot the NBR layers and NAFI masked layers

Create a colour scheme

Code 
colours <- rev(brewer.pal(11, "RdBu"))
Code 
for(i in 1:nlyr(nbr)) terra::plot(nbr[[i]], 
                                  col = colours, 
                                  # breaks = seq(-1,1, length.out = length(colours) + 1),
                                  breaks = c(-1, -0.5, -0.25, -0.1, 0.1, 0.25, 0.5, 1),
                                  main = paste0("dNBR ", time(nbr[[i]])))

Code 
for(i in 1:nlyr(nbr_NAFImask)) terra::plot(nbr_NAFImask[[i]], 
                                           col = colours, 
                                           # breaks = seq(-1,1, length.out = length(colours) + 1),
                                           breaks = c(-1, -0.5, -0.25, -0.1, 0.1, 0.25, 0.5, 1),
                                           main = paste0("dNBR NAFImask ", time(nbr_NAFImask[[i]])))

Load fitted model objects

These models were fitted with the glmmTMB package, and have a lag to the NBR covariate of 0, 1, 2 or 3 months.

Code 
dry_lag0 <- readRDS("fitted_models/buffalo_ssf10rs_dry_nbrNAFI_lag0_2024-06-13.rds")
dry_lag1 <- readRDS("fitted_models/buffalo_ssf10rs_dry_nbr_NAFImask_lag1_2024-06-13.rds")
dry_lag2 <- readRDS("fitted_models/buffalo_ssf10rs_dry_nbr_NAFImask_lag2_2024-06-13.rds")
dry_lag3 <- readRDS("fitted_models/buffalo_ssf10rs_dry_nbr_NAFImask_lag3_2024-06-13.rds")
Code 
wet_lag0 <- readRDS("fitted_models/buffalo_ssf10rs_wet_nbrNAFI_lag0_2024-06-13.rds")
wet_lag1 <- readRDS("fitted_models/buffalo_ssf10rs_wet_nbr_NAFImask_lag1_2024-06-13.rds")
wet_lag2 <- readRDS("fitted_models/buffalo_ssf10rs_wet_nbr_NAFImask_lag2_2024-06-13.rds")
wet_lag3 <- readRDS("fitted_models/buffalo_ssf10rs_wet_nbr_NAFImask_lag3_2024-06-13.rds")

Check model outputs

Dry season - NBR mask

Code 
dry_lag0 <- glmmTMB::up2date(dry_lag0)
summary(dry_lag0)
 Family: poisson  ( log )
Formula:          
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_end_scaled +  
    water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    (1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |  
    id) + (0 + nbr_NAFImask_end_scaled | id) + (0 + water_dist_end_scaled |  
    id)
Data: ssfdat_dry

      AIC       BIC    logLik  deviance  df.resid 
1158463.1 1158587.8 -579220.6 1158441.1    619424 

Random effects:

Conditional model:
 Groups  Name                    Variance  Std.Dev. 
 step_id (Intercept)             1.000e+06 1000.0000
 id      ndvi_end_scaled         4.657e-02    0.2158
 id.1    I(ndvi_end_scaled^2)    4.288e-02    0.2071
 id.2    nbr_NAFImask_end_scaled 1.101e-02    0.1049
 id.3    water_dist_end_scaled   1.944e-01    0.4409
Number of obs: 619435, groups:  step_id, 56326; id, 15

Conditional model:
                         Estimate Std. Error z value Pr(>|z|)    
ndvi_end_scaled         -0.064210   0.056860  -1.129 0.258790    
I(ndvi_end_scaled^2)    -0.201018   0.054289  -3.703 0.000213 ***
nbr_NAFImask_end_scaled  0.068109   0.031062   2.193 0.028330 *  
water_dist_end_scaled   -0.534262   0.119417  -4.474 7.68e-06 ***
sl_scaled                0.023810   0.005625   4.233 2.31e-05 ***
log_sl_scaled           -0.006310   0.005485  -1.150 0.249948    
cos_ta_scaled           -0.015570   0.004440  -3.507 0.000454 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
dry_lag1 <- glmmTMB::up2date(dry_lag1) # as model was fitted with older version of glmmTMB
summary(dry_lag1)
 Family: poisson  ( log )
Formula:          
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_lag_1_end_scaled +  
    water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    (1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |  
    id) + (0 + nbr_NAFImask_lag_1_end_scaled | id) + (0 + water_dist_end_scaled |  
    id)
Data: ssfdat_dry

      AIC       BIC    logLik  deviance  df.resid 
1158413.3 1158538.0 -579195.6 1158391.3    619324 

Random effects:

Conditional model:
 Groups  Name                          Variance  Std.Dev. 
 step_id (Intercept)                   1.000e+06 1.000e+03
 id      ndvi_end_scaled               4.862e-02 2.205e-01
 id.1    I(ndvi_end_scaled^2)          4.140e-02 2.035e-01
 id.2    nbr_NAFImask_lag_1_end_scaled 4.793e-03 6.923e-02
 id.3    water_dist_end_scaled         1.973e-01 4.442e-01
Number of obs: 619335, groups:  step_id, 56326; id, 15

Conditional model:
                               Estimate Std. Error z value Pr(>|z|)    
ndvi_end_scaled               -0.059419   0.058059  -1.023 0.306106    
I(ndvi_end_scaled^2)          -0.196315   0.053372  -3.678 0.000235 ***
nbr_NAFImask_lag_1_end_scaled  0.049287   0.021860   2.255 0.024153 *  
water_dist_end_scaled         -0.539565   0.120229  -4.488 7.20e-06 ***
sl_scaled                      0.024931   0.005623   4.434 9.26e-06 ***
log_sl_scaled                 -0.007381   0.005480  -1.347 0.178029    
cos_ta_scaled                 -0.015518   0.004440  -3.495 0.000474 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
dry_lag2 <- glmmTMB::up2date(dry_lag2) # as model was fitted with older version of glmmTMB
summary(dry_lag2)
 Family: poisson  ( log )
Formula:          
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_lag_2_end_scaled +  
    water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    (1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |  
    id) + (0 + nbr_NAFImask_lag_2_end_scaled | id) + (0 + water_dist_end_scaled |  
    id)
Data: ssfdat_dry

      AIC       BIC    logLik  deviance  df.resid 
1158327.5 1158452.2 -579152.8 1158305.5    619211 

Random effects:

Conditional model:
 Groups  Name                          Variance  Std.Dev. 
 step_id (Intercept)                   1.000e+06 1000.0000
 id      ndvi_end_scaled               5.088e-02    0.2256
 id.1    I(ndvi_end_scaled^2)          4.104e-02    0.2026
 id.2    nbr_NAFImask_lag_2_end_scaled 1.196e-02    0.1093
 id.3    water_dist_end_scaled         2.011e-01    0.4485
Number of obs: 619222, groups:  step_id, 56326; id, 15

Conditional model:
                               Estimate Std. Error z value Pr(>|z|)    
ndvi_end_scaled               -0.063999   0.059344  -1.078 0.280840    
I(ndvi_end_scaled^2)          -0.197331   0.053151  -3.713 0.000205 ***
nbr_NAFImask_lag_2_end_scaled -0.009325   0.031048  -0.300 0.763906    
water_dist_end_scaled         -0.536862   0.121287  -4.426 9.58e-06 ***
sl_scaled                      0.024958   0.005622   4.440 9.01e-06 ***
log_sl_scaled                 -0.006762   0.005483  -1.233 0.217506    
cos_ta_scaled                 -0.015806   0.004440  -3.560 0.000372 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
dry_lag3 <- glmmTMB::up2date(dry_lag3) # as model was fitted with older version of glmmTMB
summary(dry_lag3)
 Family: poisson  ( log )
Formula:          
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_lag_3_end_scaled +  
    water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    (1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |  
    id) + (0 + nbr_NAFImask_lag_3_end_scaled | id) + (0 + water_dist_end_scaled |  
    id)
Data: ssfdat_dry

      AIC       BIC    logLik  deviance  df.resid 
1154382.3 1154506.9 -577180.1 1154360.3    617135 

Random effects:

Conditional model:
 Groups  Name                          Variance  Std.Dev. 
 step_id (Intercept)                   1.000e+06 1000.0000
 id      ndvi_end_scaled               4.670e-02    0.2161
 id.1    I(ndvi_end_scaled^2)          4.085e-02    0.2021
 id.2    nbr_NAFImask_lag_3_end_scaled 2.016e-02    0.1420
 id.3    water_dist_end_scaled         1.892e-01    0.4350
Number of obs: 617146, groups:  step_id, 56273; id, 15

Conditional model:
                               Estimate Std. Error z value Pr(>|z|)    
ndvi_end_scaled               -0.067141   0.056958  -1.179 0.238489    
I(ndvi_end_scaled^2)          -0.198354   0.053026  -3.741 0.000184 ***
nbr_NAFImask_lag_3_end_scaled -0.025078   0.039013  -0.643 0.520350    
water_dist_end_scaled         -0.549929   0.117981  -4.661 3.14e-06 ***
sl_scaled                      0.025689   0.005626   4.566 4.96e-06 ***
log_sl_scaled                 -0.005391   0.005496  -0.981 0.326675    
cos_ta_scaled                 -0.015540   0.004449  -3.493 0.000478 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wet season - NBR mask

Code 
wet_lag0 <- glmmTMB::up2date(wet_lag0) # as model was fitted with older version of glmmTMB
summary(wet_lag0)
 Family: poisson  ( log )
Formula:          
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_end_scaled +  
    water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    (1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |  
    id) + (0 + nbr_NAFImask_end_scaled | id) + (0 + water_dist_end_scaled |  
    id)
Data: ssfdat_wet

      AIC       BIC    logLik  deviance  df.resid 
 780451.5  780571.9 -390214.8  780429.5    416401 

Random effects:

Conditional model:
 Groups  Name                    Variance  Std.Dev. 
 step_id (Intercept)             1.000e+06 1.000e+03
 id      ndvi_end_scaled         4.079e-02 2.020e-01
 id.1    I(ndvi_end_scaled^2)    3.577e-03 5.981e-02
 id.2    nbr_NAFImask_end_scaled 6.461e-03 8.038e-02
 id.3    water_dist_end_scaled   3.911e-01 6.254e-01
Number of obs: 416412, groups:  step_id, 39336; id, 15

Conditional model:
                         Estimate Std. Error z value Pr(>|z|)    
ndvi_end_scaled         -0.110719   0.055136  -2.008  0.04463 *  
I(ndvi_end_scaled^2)    -0.069847   0.016921  -4.128 3.66e-05 ***
nbr_NAFImask_end_scaled  0.032391   0.023856   1.358  0.17455    
water_dist_end_scaled   -0.326782   0.175689  -1.860  0.06288 .  
sl_scaled               -0.019121   0.006832  -2.799  0.00513 ** 
log_sl_scaled            0.089186   0.007292  12.231  < 2e-16 ***
cos_ta_scaled            0.029285   0.005464   5.360 8.32e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
wet_lag1 <- glmmTMB::up2date(wet_lag1) # as model was fitted with older version of glmmTMB
summary(wet_lag1)
 Family: poisson  ( log )
Formula:          
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_lag_1_end_scaled +  
    water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    (1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |  
    id) + (0 + nbr_NAFImask_lag_1_end_scaled | id) + (0 + water_dist_end_scaled |  
    id)
Data: ssfdat_wet

      AIC       BIC    logLik  deviance  df.resid 
 765469.0  765589.1 -382723.5  765447.0    408459 

Random effects:

Conditional model:
 Groups  Name                          Variance  Std.Dev. 
 step_id (Intercept)                   1.000e+06 1.000e+03
 id      ndvi_end_scaled               3.354e-02 1.831e-01
 id.1    I(ndvi_end_scaled^2)          3.271e-03 5.720e-02
 id.2    nbr_NAFImask_lag_1_end_scaled 2.537e-03 5.037e-02
 id.3    water_dist_end_scaled         3.363e-01 5.799e-01
Number of obs: 408470, groups:  step_id, 39172; id, 15

Conditional model:
                               Estimate Std. Error z value Pr(>|z|)    
ndvi_end_scaled               -0.105740   0.050571  -2.091  0.03654 *  
I(ndvi_end_scaled^2)          -0.068177   0.016291  -4.185 2.85e-05 ***
nbr_NAFImask_lag_1_end_scaled  0.083075   0.017153   4.843 1.28e-06 ***
water_dist_end_scaled         -0.205804   0.164563  -1.251  0.21108    
sl_scaled                     -0.020515   0.006906  -2.971  0.00297 ** 
log_sl_scaled                  0.090305   0.007383  12.232  < 2e-16 ***
cos_ta_scaled                  0.030390   0.005527   5.498 3.84e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
wet_lag2 <- glmmTMB::up2date(wet_lag2) # as model was fitted with older version of glmmTMB
summary(wet_lag2)
 Family: poisson  ( log )
Formula:          
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_lag_2_end_scaled +  
    water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    (1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |  
    id) + (0 + nbr_NAFImask_lag_2_end_scaled | id) + (0 + water_dist_end_scaled |  
    id)
Data: ssfdat_wet

      AIC       BIC    logLik  deviance  df.resid 
 767710.8  767831.0 -383844.4  767688.8    409719 

Random effects:

Conditional model:
 Groups  Name                          Variance  Std.Dev. 
 step_id (Intercept)                   1.000e+06 1.000e+03
 id      ndvi_end_scaled               3.798e-02 1.949e-01
 id.1    I(ndvi_end_scaled^2)          3.488e-03 5.906e-02
 id.2    nbr_NAFImask_lag_2_end_scaled 1.282e-02 1.132e-01
 id.3    water_dist_end_scaled         2.998e-01 5.475e-01
Number of obs: 409730, groups:  step_id, 39270; id, 15

Conditional model:
                               Estimate Std. Error z value Pr(>|z|)    
ndvi_end_scaled               -0.103767   0.053428  -1.942 0.052115 .  
I(ndvi_end_scaled^2)          -0.066154   0.016748  -3.950 7.81e-05 ***
nbr_NAFImask_lag_2_end_scaled  0.121523   0.032968   3.686 0.000228 ***
water_dist_end_scaled         -0.177165   0.156458  -1.132 0.257485    
sl_scaled                     -0.017556   0.006867  -2.557 0.010570 *  
log_sl_scaled                  0.094926   0.007394  12.839  < 2e-16 ***
cos_ta_scaled                  0.032153   0.005519   5.825 5.70e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code 
wet_lag3 <- glmmTMB::up2date(wet_lag3) # as model was fitted with older version of glmmTMB
summary(wet_lag3)
 Family: poisson  ( log )
Formula:          
case_ ~ -1 + ndvi_end_scaled + I(ndvi_end_scaled^2) + nbr_NAFImask_lag_3_end_scaled +  
    water_dist_end_scaled + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    (1 | step_id) + (0 + ndvi_end_scaled | id) + (0 + I(ndvi_end_scaled^2) |  
    id) + (0 + nbr_NAFImask_lag_3_end_scaled | id) + (0 + water_dist_end_scaled |  
    id)
Data: ssfdat_wet

      AIC       BIC    logLik  deviance  df.resid 
 767616.8  767736.9 -383797.4  767594.8    409701 

Random effects:

Conditional model:
 Groups  Name                          Variance  Std.Dev. 
 step_id (Intercept)                   1.000e+06 1.000e+03
 id      ndvi_end_scaled               3.876e-02 1.969e-01
 id.1    I(ndvi_end_scaled^2)          3.630e-03 6.025e-02
 id.2    nbr_NAFImask_lag_3_end_scaled 1.429e-03 3.780e-02
 id.3    water_dist_end_scaled         2.790e-01 5.282e-01
Number of obs: 409712, groups:  step_id, 39298; id, 15

Conditional model:
                               Estimate Std. Error z value Pr(>|z|)    
ndvi_end_scaled               -0.107474   0.053884  -1.995 0.046092 *  
I(ndvi_end_scaled^2)          -0.065768   0.017046  -3.858 0.000114 ***
nbr_NAFImask_lag_3_end_scaled  0.063826   0.016883   3.780 0.000157 ***
water_dist_end_scaled         -0.199873   0.151718  -1.317 0.187705    
sl_scaled                     -0.017101   0.006864  -2.491 0.012723 *  
log_sl_scaled                  0.092028   0.007384  12.463  < 2e-16 ***
cos_ta_scaled                  0.031568   0.005519   5.720 1.07e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Create dataframes to plot the coefficients

Code 
coef_df_dry_lag0 <- data.frame("model" = "dry_lag0", 
                      "model_covariate" = names(fixef(dry_lag0)$cond),
                      "Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
                      "Estimate" = coef(summary(dry_lag0))$cond[, "Estimate"],
                      "SE" = coef(summary(dry_lag0))$cond[, "Std. Error"]
                      ) 

coef_df_dry_lag1 <- data.frame("model" = "dry_lag1", 
                      "model_covariate" = names(fixef(dry_lag1)$cond),
                      "Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
                      "Estimate" = coef(summary(dry_lag1))$cond[, "Estimate"],
                      "SE" = coef(summary(dry_lag1))$cond[, "Std. Error"]
                      )

coef_df_dry_lag2 <- data.frame("model" = "dry_lag2", 
                      "model_covariate" = names(fixef(dry_lag2)$cond),
                      "Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
                      "Estimate" = coef(summary(dry_lag2))$cond[, "Estimate"],
                      "SE" = coef(summary(dry_lag2))$cond[, "Std. Error"]
                      )

coef_df_dry_lag3 <- data.frame("model" = "dry_lag3", 
                      "model_covariate" = names(fixef(dry_lag3)$cond),
                      "Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
                      "Estimate" = coef(summary(dry_lag3))$cond[, "Estimate"],
                      "SE" = coef(summary(dry_lag3))$cond[, "Std. Error"]
                      )

coef_df_dry <- rbind(coef_df_dry_lag0, coef_df_dry_lag1, coef_df_dry_lag2, coef_df_dry_lag3)

coef_df_dry <- coef_df_dry %>%
  
  mutate(
    # Calculate confidence intervals
    LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE, 
    UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
    LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
    UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
    LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
    UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
    # Compute p-values
    pvalue = 2 * pnorm(-abs(Estimate) / SE)
  )

# Add stars indicating the significance
coef_df_dry$Significance <- sapply(1:nrow(coef_df_dry), function(x){
  if (coef_df_dry$pvalue[x] <= 0.001){
    return("***")
  } else if (coef_df_dry$pvalue[x] <= 0.01) {
    return("**")
  } else if (coef_df_dry$pvalue[x] <= 0.05) {
    return("*")
  }
})

# Remove the intercept term
coef_df_dry <- coef_df_dry %>% filter(Covariate != "Intercept")

# Add a column indicating the preference
coef_df_dry$Preference <- ifelse(coef_df_dry$Estimate > 0, "Preferred", "Avoided")
coef_df_dry$Preference <- factor(coef_df_dry$Preference, levels = c("Avoided", "Preferred"))

Prepare data frame for plotting

Specify the order in which the coefficients should be plotted

Code 
order <- c(
  "NDVI",
  "NDVI^2",
  "NBR",
  "DistanceToWater",
  "sl",
  "log(sl)",
  "cos(ta)"
)

Prepare dataset for plotting confidence intervals

Code 
# Reverse the order of the 'model' factor levels
coef_df_dry$model <- factor(coef_df_dry$model, levels = rev(levels(factor(coef_df_dry$model))))

coef_df_dry2 <- coef_df_dry %>%
  dplyr::select(model, Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

Prepare plot with covariates on the y-axis and the estimates on the x-axis

Code 
dry_nbr_lag_plot <- ggplot(data = coef_df_dry, 
                           aes(y = Covariate, x = Estimate, col = factor(Preference), group = factor(model))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  annotate(geom = "segment",
           x      = -0.85, xend = 0.25,
           y      = 3.5, yend   = 3.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5, position = position_dodge(width = 0.5)) +
  geom_errorbarh(data = coef_df_dry2,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5, position = position_dodge(width = 0.5)) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F, position = position_dodge(width = 0.5)) +
  scale_y_discrete(limits = rev(order)) +
  theme_classic() +
  scale_x_continuous(limits = c(-0.9, 0.25)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#FF8123", "#9B4200")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#FF8123")))

dry_nbr_lag_plot
Warning: Removed 10 rows containing missing values or values outside the scale range
(`geom_text()`).

Code 
# ggsave("outputs/plots/nbr_lag_model_dry.png", device = "png", 
#        width = 150, height = 110, units = "mm", scale = 1, dpi = 1000)

Most ‘informative’ model

Code 
dry_nbr_lag_plot <- ggplot(data = coef_df_dry %>% filter(model == "dry_lag0"), 
                           aes(y = Covariate, x = Estimate, col = factor(Preference), group = factor(model))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  annotate(geom = "segment",
           x      = -0.85, xend = 0.25,
           y      = 3.5, yend   = 3.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5, position = position_dodge(width = 0.5)) +
  geom_errorbarh(data = coef_df_dry2 %>% filter(model == "dry_lag0"),
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5, position = position_dodge(width = 0.5)) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F, position = position_dodge(width = 0.5)) +
  scale_y_discrete(limits = rev(order)) +
  theme_classic() +
  scale_x_continuous(limits = c(-0.9, 0.25)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#FF8123", "#9B4200")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#FF8123")))

dry_nbr_lag_plot
Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_text()`).

Code 
# ggsave("outputs/plots/nbr_lag_model_final_dry.png", device = "png", 
#        width = 150, height = 70, units = "mm", scale = 1, dpi = 1000)

NDVI quadratic curve

Pull out the coefficients from the fitted model.

Code 
best_nbr_model <- coef_df_dry %>% filter(model == "dry_lag0")
ndvi_param <- best_nbr_model %>% filter(Covariate == "NDVI") %>% pull(Estimate) / ndvi_scale_sd
ndvi_param_se <- best_nbr_model %>% filter(Covariate == "NDVI") %>% pull(SE) / ndvi_scale_sd
ndvi2_param <- best_nbr_model %>% filter(Covariate == "NDVI^2") %>% pull(Estimate) / ndvi_scale_sd
ndvi2_param_se <- best_nbr_model %>% filter(Covariate == "NDVI^2") %>% pull(SE) / ndvi_scale_sd

Create a dataframe of x and y values for the quadratic curve.

Code 
ndvi_df <- data.frame("NDVI" = seq(-0.2, 0.8, 0.01))
ndvi_df$NDVI2 <- ndvi_df$NDVI^2
ndvi_df$response <- (ndvi_df$NDVI * ndvi_param) + (ndvi_df$NDVI2 * ndvi2_param)
ndvi_df$RSS <- exp(ndvi_df$response)

Plot the quadratic curve

Code 
ggplot() +
  geom_hline(yintercept = 1, color = "black", lty = "dashed", lwd = 0.3) +
  geom_line(data = ndvi_df, aes(x = NDVI, y = RSS), color = "#FF8123", size = 1) +
  scale_y_continuous("Relative Selection Strength (RSS)") +
  # scale_x_continuous("NDVI", limits = c(-0.1, 0.75)) +
  theme_classic() 
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Code 
# ggsave("outputs/plots/ndvi_quadratic_curve_dry.png", device = "png", 
#        width = 150, height = 70, units = "mm", scale = 1, dpi = 1000)

Calculate confidence intervals for the curve.

Here we simulate the beta values for NDVI and NDVI^2, and calculate the RSS for each combination of beta values and NDVI values. We then calculate the 95% confidence intervals for the RSS values.

Code 
# Number of simulations
n_sim <- 10000

# Simulate beta values
ndvi_sim <- rnorm(n_sim, mean = ndvi_param, sd = ndvi_param_se)
# hist(ndvi_sim)
ndvi2_sim <- rnorm(n_sim, mean = ndvi2_param, sd = ndvi2_param_se)
# hist(ndvi2_sim)

# Define x values
ndvi_values <- seq(-0.2, 0.8, length.out = 100)

# Initialize matrix to store RSS simulations
RSS_sim_matrix <- matrix(NA, nrow = n_sim, ncol = length(ndvi_values))

# Calculate RSS for each simulated beta and x
for (i in 1:length(ndvi_values)) {
    RSS_sim_matrix[, i] <- exp(ndvi_sim * ndvi_values[i] + ndvi2_sim * ndvi_values[i]^2)
}

# Calculate confidence intervals for each x
RSS_CI_lower <- apply(RSS_sim_matrix, 2, quantile, probs = 0.025)
RSS_CI_upper <- apply(RSS_sim_matrix, 2, quantile, probs = 0.975)

# Convert the matrix to a data frame
RSS_sim_df <- as.data.frame(RSS_sim_matrix)

# Add column names corresponding to x_values
colnames(RSS_sim_df) <- paste0("x_", ndvi_values)

# Add a simulation ID column
RSS_sim_df$Simulation <- 1:n_sim

# Reshape the data frame from wide to long format
RSS_long_df <- pivot_longer(
  RSS_sim_df,
  cols = starts_with("x_"),
  names_to = "x",
  values_to = "RSS_sim"
)

# Remove the "x_" prefix and convert x to numeric
RSS_long_df$x <- as.numeric(sub("x_", "", RSS_long_df$x))

# Calculate mean and confidence intervals for each x
RSS_summary <- RSS_long_df %>%
  group_by(x) %>%
  summarise(
    RSS_mean = mean(RSS_sim),
    Lower_CI = quantile(RSS_sim, probs = 0.025),
    Lower_CI_50 = quantile(RSS_sim, probs = 0.25),
    Upper_CI = quantile(RSS_sim, probs = 0.975),
    Upper_CI_50 = quantile(RSS_sim, probs = 0.75)
  )

# Plot mean RSS with confidence intervals
ggplot(RSS_summary, aes(x = x, y = RSS_mean)) +
  geom_hline(yintercept = 1, color = "black", lty = "dashed", lwd = 0.3) +
  geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), fill = "#FF8123", alpha = 0.15) +
  geom_ribbon(aes(ymin = Lower_CI_50, ymax = Upper_CI_50), fill = "#FF8123", alpha = 0.15) +
  geom_line(color = "#FF8123") +
  labs(
    # title = "Relative Selection Strength with 50% and 95% Confidence Intervals",
    x = "NDVI",
    y = "Relative Selection Strength (RSS)"
  ) +
  theme_bw()

Code 
ggsave("outputs/plots/ndvi_quadratic_curve_CI_dry.png", device = "png",
       width = 150, height = 100, units = "mm", scale = 1, dpi = 1000)
Code 
coef_df_wet_lag0 <- data.frame("model" = "wet_lag0", 
                      "model_covariate" = names(fixef(wet_lag0)$cond),
                      "Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
                      "Estimate" = coef(summary(wet_lag0))$cond[, "Estimate"],
                      "SE" = coef(summary(wet_lag0))$cond[, "Std. Error"]
                      ) 

coef_df_wet_lag1 <- data.frame("model" = "wet_lag1", 
                      "model_covariate" = names(fixef(wet_lag1)$cond),
                      "Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
                      "Estimate" = coef(summary(wet_lag1))$cond[, "Estimate"],
                      "SE" = coef(summary(wet_lag1))$cond[, "Std. Error"]
                      )

coef_df_wet_lag2 <- data.frame("model" = "wet_lag2", 
                      "model_covariate" = names(fixef(wet_lag2)$cond),
                      "Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
                      "Estimate" = coef(summary(wet_lag2))$cond[, "Estimate"],
                      "SE" = coef(summary(wet_lag2))$cond[, "Std. Error"]
                      )

coef_df_wet_lag3 <- data.frame("model" = "wet_lag3", 
                      "model_covariate" = names(fixef(wet_lag3)$cond),
                      "Covariate" = c("NDVI", "NDVI^2", "NBR", "DistanceToWater", "sl", "log(sl)", "cos(ta)"),
                      "Estimate" = coef(summary(wet_lag3))$cond[, "Estimate"],
                      "SE" = coef(summary(wet_lag3))$cond[, "Std. Error"]
                      )

coef_df_wet <- rbind(coef_df_wet_lag0, coef_df_wet_lag1, coef_df_wet_lag2, coef_df_wet_lag3)

coef_df_wet <- coef_df_wet %>%
  
  mutate(
    # Calculate confidence intervals
    LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE, 
    UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
    LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
    UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
    LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
    UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
    # Compute p-values
    pvalue = 2 * pnorm(-abs(Estimate) / SE)
  )

# Add stars indicating the significance
coef_df_wet$Significance <- sapply(1:nrow(coef_df_wet), function(x){
  if (coef_df_wet$pvalue[x] <= 0.001){
    return("***")
  } else if (coef_df_wet$pvalue[x] <= 0.01) {
    return("**")
  } else if (coef_df_wet$pvalue[x] <= 0.05) {
    return("*")
  }
})

# Remove the intercept term
coef_df_wet <- coef_df_wet %>% filter(Covariate != "Intercept")

# Add a column indicating the preference
coef_df_wet$Preference <- ifelse(coef_df_wet$Estimate > 0, "Preferred", "Avoided")
coef_df_wet$Preference <- factor(coef_df_wet$Preference, levels = c("Avoided", "Preferred"))

Prepare confidence interval dataframe

Code 
# Reverse the order of the 'model' factor levels
coef_df_wet$model <- factor(coef_df_wet$model, levels = rev(levels(factor(coef_df_wet$model))))

coef_df_wet2 <- coef_df_wet %>%
  dplyr::select(model, Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

Prepare plot with covariates on the y-axis and the estimates on the x-axis

Code 
wet_nbr_lag_plot <- ggplot(data = coef_df_wet, 
                           aes(y = Covariate, x = Estimate, col = factor(Preference), group = factor(model))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  annotate(geom = "segment",
           x      = -0.9, xend = 0.25,
           y      = 3.5, yend   = 3.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5, position = position_dodge(width = 0.5)) +
  geom_errorbarh(data = coef_df_wet2,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5, position = position_dodge(width = 0.5)) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F, position = position_dodge(width = 0.5)) +
  scale_y_discrete(limits = rev(order)) +
  theme_classic() +
  scale_x_continuous(limits = c(-0.9, 0.25)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#238DFF", "#004691")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#1377E2")))

wet_nbr_lag_plot
Warning: Removed 6 rows containing missing values or values outside the scale range
(`geom_text()`).

Code 
# ggsave("outputs/plots/nbr_lag_model_wet.png", device = "png", 
#        width = 150, height = 110, units = "mm", scale = 1, dpi = 1000)

Most ‘informative’ model

Code 
wet_nbr_lag_plot <- ggplot(data = coef_df_wet %>% filter(model == "wet_lag2"), 
                           aes(y = Covariate, x = Estimate, col = factor(Preference), group = factor(model))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  annotate(geom = "segment",
           x      = -0.9, xend = 0.25,
           y      = 3.5, yend   = 3.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5, position = position_dodge(width = 0.5)) +
  geom_errorbarh(data = coef_df_wet2 %>% filter(model == "wet_lag2"),
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5, position = position_dodge(width = 0.5)) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F, position = position_dodge(width = 0.5)) +
  scale_y_discrete(limits = rev(order)) +
  theme_classic() +
  scale_x_continuous(limits = c(-0.9, 0.25)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#238DFF", "#004691")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#1377E2")))

wet_nbr_lag_plot
Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_text()`).

Code 
# ggsave("outputs/plots/nbr_lag_model_final_wet.png", device = "png", 
#        width = 150, height = 70, units = "mm", scale = 1, dpi = 1000)

NDVI quadratic curve

Pull out the coefficients from the fitted model.

Code 
best_nbr_model <- coef_df_wet %>% filter(model == "wet_lag0")
ndvi_param <- best_nbr_model %>% filter(Covariate == "NDVI") %>% pull(Estimate) / ndvi_scale_sd
ndvi_param_se <- best_nbr_model %>% filter(Covariate == "NDVI") %>% pull(SE) / ndvi_scale_sd
ndvi2_param <- best_nbr_model %>% filter(Covariate == "NDVI^2") %>% pull(Estimate) / ndvi_scale_sd
ndvi2_param_se <- best_nbr_model %>% filter(Covariate == "NDVI^2") %>% pull(SE) / ndvi_scale_sd

Create a dataframe of x and y values for the quadratic curve.

Code 
ndvi_df <- data.frame("NDVI" = seq(-0.2, 0.8, 0.01))
ndvi_df$NDVI2 <- ndvi_df$NDVI^2
ndvi_df$response <- (ndvi_df$NDVI * ndvi_param) + (ndvi_df$NDVI2 * ndvi2_param)
ndvi_df$RSS <- exp(ndvi_df$response)

Plot the quadratic curve.

Code 
ggplot() +
  geom_hline(yintercept = 1, color = "black", lty = "dashed", lwd = 0.3) +
  geom_line(data = ndvi_df, aes(x = NDVI, y = RSS), color = "#1377E2", size = 1) +
  scale_y_continuous("Relative Selection Strength (RSS)") +
  # scale_x_continuous("NDVI", limits = c(-0.25, 0.75)) +
  theme_classic()

Code 
# ggsave("outputs/plots/ndvi_quadratic_curve_wet.png", device = "png", 
#        width = 150, height = 70, units = "mm", scale = 1, dpi = 1000)

Calculate confidence intervals for the curve.

Here we simulate the beta values for NDVI and NDVI^2, and calculate the RSS for each combination of beta values and NDVI values. We then calculate the 95% confidence intervals for the RSS values.

Code 
# Number of simulations
n_sim <- 10000

# Simulate beta values
ndvi_sim <- rnorm(n_sim, mean = ndvi_param, sd = ndvi_param_se)
# hist(ndvi_sim)
ndvi2_sim <- rnorm(n_sim, mean = ndvi2_param, sd = ndvi2_param_se)
# hist(ndvi2_sim)

# Define x values
ndvi_values <- seq(-0.2, 0.8, length.out = 100)

# Initialize matrix to store RSS simulations
RSS_sim_matrix <- matrix(NA, nrow = n_sim, ncol = length(ndvi_values))

# Calculate RSS for each simulated beta and x
for (i in 1:length(ndvi_values)) {
    RSS_sim_matrix[, i] <- exp(ndvi_sim * ndvi_values[i] + ndvi2_sim * ndvi_values[i]^2)
}

# Calculate confidence intervals for each x
RSS_CI_lower <- apply(RSS_sim_matrix, 2, quantile, probs = 0.025)
RSS_CI_upper <- apply(RSS_sim_matrix, 2, quantile, probs = 0.975)

# Convert the matrix to a data frame
RSS_sim_df <- as.data.frame(RSS_sim_matrix)

# Add column names corresponding to x_values
colnames(RSS_sim_df) <- paste0("x_", ndvi_values)

# Add a simulation ID column
RSS_sim_df$Simulation <- 1:n_sim

# Reshape the data frame from wide to long format
RSS_long_df <- pivot_longer(
  RSS_sim_df,
  cols = starts_with("x_"),
  names_to = "x",
  values_to = "RSS_sim"
)

# Remove the "x_" prefix and convert x to numeric
RSS_long_df$x <- as.numeric(sub("x_", "", RSS_long_df$x))

# Calculate mean and confidence intervals for each x
RSS_summary <- RSS_long_df %>%
  group_by(x) %>%
  summarise(
    RSS_mean = mean(RSS_sim),
    Lower_CI = quantile(RSS_sim, probs = 0.025),
    Lower_CI_50 = quantile(RSS_sim, probs = 0.25),
    Upper_CI = quantile(RSS_sim, probs = 0.975),
    Upper_CI_50 = quantile(RSS_sim, probs = 0.75)
  )

# Plot mean RSS with confidence intervals
ggplot(RSS_summary, aes(x = x, y = RSS_mean)) +
  geom_hline(yintercept = 1, color = "black", lty = "dashed", lwd = 0.3) +
  geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), fill = "#1377E2", alpha = 0.15) +
  geom_ribbon(aes(ymin = Lower_CI_50, ymax = Upper_CI_50), fill = "#1377E2", alpha = 0.15) +
  geom_line(color = "#1377E2") +
  labs(
    # title = "Relative Selection Strength with 50% and 95% Confidence Intervals",
    x = "NDVI",
    y = "Relative Selection Strength (RSS)"
  ) +
  theme_bw()

Code 
ggsave("outputs/plots/ndvi_quadratic_curve_CI_wet.png", device = "png",
       width = 150, height = 100, units = "mm", scale = 1, dpi = 1000)

Checking used vs avilable distributions

NDVI

Code 
# bin the ssf_data$nbr_end values
# ssf_data$ndvi_end_bin <- cut(ssf_data$ndvi_end, breaks = c(-1, -0.5, -0.25, -0.1, 0.1, 0.25, 0.5, 1))
ssf_data$ndvi_end_bin <- cut(ssf_data$ndvi_end, breaks = seq(-1,1,0.1))
ssf_data_plot <- ssf_data %>% filter(ndvi_end > -0.2 & ndvi_end < 0.8)

# inspect the data for dry season 
ssf_data_plot %>% 
  filter(season == "dry") %>%
  drop_na() %>% 
  group_by(case_, ndvi_end_bin) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(ndvi_end_bin, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "NDVI", y = "Proportion", fill = "Case")+
  ggtitle("Dry season")+
  scale_fill_brewer(palette = "Paired", name="case_", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Code 
# inspect the data for wet season 
ssf_data_plot %>% 
  filter(season == "wet") %>%
  drop_na() %>% 
  group_by(case_, ndvi_end_bin) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(ndvi_end_bin, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "NDVI", y = "Proportion", fill = "Case")+
  ggtitle("Wet season")+
  scale_fill_brewer(palette = "Paired", name="case_", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

NBR

Code 
# bin the ssf_data$nbr_end values
ssf_data$nbr_lag_0_end_bin <- cut(ssf_data$nbr_lag_0_end, breaks = c(-1, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 1))

# inspect the data for dry season 
ssf_data %>% 
  filter(season == "dry") %>%
  drop_na() %>% 
  group_by(case_, nbr_lag_0_end_bin) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(nbr_lag_0_end_bin, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "Case")+
  ggtitle("Dry season - no NBR lag")+
  scale_fill_brewer(palette = "Paired", name="case_", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Code 
# inspect the data for wet season 
ssf_data %>% 
  filter(season == "wet") %>%
  drop_na() %>% 
  group_by(case_, nbr_lag_0_end_bin) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(nbr_lag_0_end_bin, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "Case")+
  ggtitle("Wet season - no NBR lag")+
  scale_fill_brewer(palette = "Paired", name="case_", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Code 
# bin the ssf_data$nbr_end values
ssf_data$nbr_lag_1_end_bin <- cut(ssf_data$nbr_lag_1_end, breaks = c(-1, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 1))

# inspect the data for dry season 
ssf_data %>% 
  filter(season == "dry") %>%
  drop_na() %>% 
  group_by(case_, nbr_lag_1_end_bin) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(nbr_lag_1_end_bin, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "Case")+
  ggtitle("Dry season - 1 month NBR lag")+
  scale_fill_brewer(palette = "Paired", name="case_", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Code 
# inspect the data for wet season 
ssf_data %>% 
  filter(season == "wet") %>%
  drop_na() %>% 
  group_by(case_, nbr_lag_1_end_bin) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(nbr_lag_1_end_bin, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "Case")+
  ggtitle("Wet season - 1 month NBR lag")+
  scale_fill_brewer(palette = "Paired", name="case_", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Code 
# bin the ssf_data$nbr_end values
ssf_data$nbr_lag_2_end_bin <- cut(ssf_data$nbr_lag_2_end, breaks = c(-1, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 1))

# inspect the data for dry season 
ssf_data %>% 
  filter(season == "dry") %>%
  drop_na() %>% 
  group_by(case_, nbr_lag_2_end_bin) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(nbr_lag_2_end_bin, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "Case")+
  ggtitle("Dry season - 2 month NBR lag")+
  scale_fill_brewer(palette = "Paired", name="case_", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Code 
# inspect the data for wet season 
ssf_data %>% 
  filter(season == "wet") %>%
  drop_na() %>% 
  group_by(case_, nbr_lag_2_end_bin) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(nbr_lag_2_end_bin, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "Case")+
  ggtitle("Wet season - 2 month NBR lag")+
  scale_fill_brewer(palette = "Paired", name="case_", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Code 
# bin the ssf_data$nbr_end values
ssf_data$nbr_lag_3_end_bin <- cut(ssf_data$nbr_lag_3_end, breaks = c(-1, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 1))

# inspect the data for dry season 
ssf_data %>% 
  filter(season == "dry") %>%
  drop_na() %>% 
  group_by(case_, nbr_lag_3_end_bin) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(nbr_lag_3_end_bin, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "Case")+
  ggtitle("Dry season - 3 month NBR lag")+
  scale_fill_brewer(palette = "Paired", name="case_", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Code 
# inspect the data for wet season 
ssf_data %>% 
  filter(season == "wet") %>%
  drop_na() %>% 
  group_by(case_, nbr_lag_3_end_bin) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(nbr_lag_3_end_bin, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "Case")+
  ggtitle("Wet season - 3 month NBR lag")+
  scale_fill_brewer(palette = "Paired", name="case_", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Code 
sessionInfo()
R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=English_Australia.utf8  LC_CTYPE=English_Australia.utf8   
[3] LC_MONETARY=English_Australia.utf8 LC_NUMERIC=C                      
[5] LC_TIME=English_Australia.utf8    

time zone: Australia/Brisbane
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] lemon_0.4.9        RColorBrewer_1.1-3 terra_1.7-78       glmmTMB_1.1.10    
 [5] lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4       
 [9] purrr_1.0.2        readr_2.1.5        tidyr_1.3.1        tibble_3.2.1      
[13] ggplot2_3.5.1      tidyverse_2.0.0   

loaded via a namespace (and not attached):
 [1] gtable_0.3.5        TMB_1.9.15          xfun_0.47          
 [4] htmlwidgets_1.6.4   lattice_0.22-6      tzdb_0.4.0         
 [7] numDeriv_2016.8-1.1 vctrs_0.6.5         tools_4.4.1        
[10] Rdpack_2.6.1        generics_0.1.3      parallel_4.4.1     
[13] fansi_1.0.6         pkgconfig_2.0.3     Matrix_1.7-0       
[16] lifecycle_1.0.4     farver_2.1.2        compiler_4.4.1     
[19] textshaping_0.4.0   munsell_0.5.1       codetools_0.2-20   
[22] htmltools_0.5.8.1   yaml_2.3.10         crayon_1.5.3       
[25] pillar_1.9.0        nloptr_2.1.1        MASS_7.3-64        
[28] reformulas_0.3.0    boot_1.3-30         nlme_3.1-164       
[31] tidyselect_1.2.1    digest_0.6.37       mvtnorm_1.3-3      
[34] stringi_1.8.4       labeling_0.4.3      splines_4.4.1      
[37] fastmap_1.2.0       grid_4.4.1          colorspace_2.1-1   
[40] cli_3.6.3           magrittr_2.0.3      utf8_1.2.4         
[43] withr_3.0.1         scales_1.3.0        bit64_4.0.5        
[46] timechange_0.3.0    estimability_1.5.1  rmarkdown_2.28     
[49] emmeans_1.10.7      bit_4.0.5           lme4_1.1-35.5      
[52] gridExtra_2.3       ragg_1.3.3          hms_1.1.3          
[55] coda_0.19-4.1       evaluate_1.0.0      knitr_1.48         
[58] rbibutils_2.2.16    mgcv_1.9-1          rlang_1.1.4        
[61] Rcpp_1.0.13         xtable_1.8-4        glue_1.7.0         
[64] vroom_1.6.5         rstudioapi_0.16.0   minqa_1.2.8        
[67] jsonlite_1.8.8      plyr_1.8.9          R6_2.5.1           
[70] systemfonts_1.1.0